Estimation of solitary pulmonary nodule diameters with reaction-diffusion segmentation

ABSTRACT

A reactive-diffusion method for estimating a diameter of an object of interest includes providing a volume of interest including a plurality of voxels, initializing at least two volumes of the volume of interest, wherein each of the voxels has at least two values corresponding to the at least two volumes of the volume of interest, respectively, performing a diffusion operation and a reaction operation on the voxels to adjust the at least two values, comparing, for each voxel, the at least two values to a threshold to assign each voxel to one of the at least two volumes, wherein the assignment of the voxels is a segmentation result, and estimating a diameter of the object of interest from the segmentation result, wherein the object of interest is represented by at least one of the at least two volumes but less than all the volumes.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of Provisional Application No.60/669,415 filed on Apr. 8, 2005 in the United States Patent andTrademark Office, the contents of which are herein incorporated byreference in its entirety.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present disclosure relates to image processing, and moreparticularly to a system and method for estimating a size of solitaryobjects in of interest in medical images.

2. Description of Related Art

Measuring the size of pulmonary nodules from X-ray computed tomography(CT) data is an important practice for diagnosis and progressionanalysis of lung cancer. The nodule size often plays an important rolein choosing a proper patient care, and is also an effective feature toseparate true nodules form nodule-like spurious findings. Typically, thesize is represented by the diameter of the nodule. Automating this taskfor computer-aided diagnosis (CAD) is, however, a difficult problem dueto intensity variations, partial volume effects, attachment to otherstructures, and noises.

A Ct-based screening protocol specified by the International Early LungCancer Action Program (I-ELCAP) details how the diameter of pulmonarynodules should be measured and how the measurements should be used fordetermining the patient management. According to the protocol, theresult of an initial CT screening of lung is considered positive if atleast one solid or part-solid nodule with 5.00 mm or more in diameter orat least one non-solid nodule with 8.0 mm or more in diameter is found.Although these 5 mm and 8 mm thresholds are likely to drop as moreaccurate screening becomes possible with high resolution multi-detectorhelical CT (MDCT), the importance of module size in cancer diagnosiswill stay unchanged.

To this end, an automated size estimation algorithm using a GaussianEllipsoid Fit (EF) has been developed. Given a marker positioned near anodule, the algorithm computes the location, orientation and radii of anellipsoid that models closely the intensity variation nearby the marker.It employs mean-shift and scale estimator to find the solution. Thevolume and diameter of the nodule can be estimated from the ellipsoid.EF can be incorporated into a CAD system where markers and provided fromeither manual markings of a reader or the detection algorithm of thesystem. The accuracy of the estimates has been verified using a largedatabase with manual size measurements. However, the technique tends tohave difficulties with small nodules due mainly to the small sample sizeproblem.

SUMMARY OF THE INVENTION

According to an embodiment of the present disclosure, areactive-diffusion method for estimating a diameter of an object ofinterest comprises providing a volume of interest comprising a pluralityof voxels and including the object of interest, and initializing atleast two volumes of the volume of interest, wherein each of theplurality of voxels has at least two values corresponding to the atleast two volumes of the volume of interest, respectively. The methodfurther includes performing, iteratively, a diffusion operation and areaction operation on the plurality of voxels to adjust the at least twovalues, comparing, for each voxel, the at least two values to athreshold to assign each voxel to one of the at least two volumes,wherein the assignment of the plurality of voxels is a segmentationresult, and estimating a diameter of the object of interest from thesegmentation result, wherein the object of interest is represented by atleast one of the at least two volumes but less than all the volumes.

BRIEF DESCRIPTION OF THE DRAWINGS

Preferred embodiments of the present invention will be described belowin more detail, with reference to the accompanying drawings:

FIG. 1A is a flow chart of a method for reactive-diffusion segmentationaccording to an embodiment of the present disclosure;

FIG. 1B is a flow chart of a method for hybrid segmentation according toan embodiment of the present disclosure;

FIGS. 2A-B are graphs of an estimate of performance of EF, LDM, and RDfor solitary nodules according to an embodiment of the presentdisclosure;

FIG. 3A-B are graphs of an estimate of performance of EF, LDM, and RDfor non-solitary nodules according to an embodiment of the presentdisclosure;

FIG. 4A-B are graphs of an estimate of performance of EF, LDM, and RDfor solitary and non-solitary nodules according to an embodiment of thepresent disclosure;

FIG. 5 is a diagram of a computer system for implementing a methodaccording to an embodiment of the present disclosure.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

According to an embodiment of the present disclosure, a method forextracting a pulmonary nodule and estimating its diameter from helicalthoracic CT scans. A rough location of the nodule is assumed to be givenby either a reader or a computer aided diagnosis program. A combinationof two segmentation methods: reaction-diffusion system (RD) andellipsoid fit (EF) may be implemented. RD is used on solitary nodulesand EF of non-solitary nodules. The solitary/non-solitary nodule type isdetermined with manual measurements of over 1300 nodules taken form over240 CT volumes. The performance of the hybrid approach is compared witha local density maximum algorithm, RD, and EF. The experiments show thatthe hybrid technique provides the most accurate size estimates, andreduces the computation time of EF by 50%.

Segmentation

Referring to FIG. 1, according to an embodiment of the presentdisclosure, a reactive-diffusion system implements a segmentationalgorithm that labels each voxel of an inputted volume of interest (VOI)(see block 101) to either foreground or background. Although a systemmay implemented more than two labels, the discussion is of a two-labelsystem for simplicity. With RD segmentation, each voxel of the inputvolume to the segmentation accompanies two values: x₀(s) and x₁(s) wheres denotes the location of the voxel. They take continuous values withthe following constraints: x₀(s)≧0,x₁(s)≧0, and x₀(s)+x₁(s)=1 for all s.These are background and foreground labels, respectively. A set ofbackground labels forms a background volume, and a set of foregroundlabels forms a foreground volume.

The background and foreground volumes are initialized at block 102. Theinitialization scheme is task dependent, and for the nodule segmentationtask, we choose the following formula for each label.x ₁(1)=e ^(−2(1−min(1,I(s)/1000))) ²   (1)and x₀(s)=1−x₁(s) where I(s) denotes the CT value of the input volume atthe locations s, and 1000 is near the upper CT value of non-calcifiedpulmonary nodules. Hence for voxels with I≧1000, x₁=1 and x₀=0. Forvoxels with I<1000, x₁ decreases monotonically as I deviates from 1000while x₀ increases. x₁ takes the minimum value of e⁻² when I=0.

After the initialization 102, the background and foreground volumesundergo diffusion 103 and reaction 104 phases alternately for a numberof times. During the diffusion phase, these volumes are diffusedseparately over the spatial domain. A linear diffusion 103 may be used,such as,x _(i) Ý(s)=λ∇² x _(i)(s)  (2)where iε{0,1} and the diffusion rate λ is set to 0.5.

During the reaction phase 104, x₀(s) and x₁(s) compete against eachother. The competition is implemented in a form of replicator dynamicsas shown below.x _(i) Ý(s)=μx _(i)(s)((Ax(s))i−x ^(T)(s)Ax(s)))   (3)where iε{0,1}, x(s)=(x₀(s)x₁(s))^(T), μ is the reaction rate,(Ax(s))_(i) denotes ith row of Ax(s), and A is a fitness matrix. Forexperiments, a 2×2 identity matrix was used as A. More complex matricesmay be used for other tasks. The update rate is set to x(s)^(T)Ax(s))⁻¹.Then,x _(i)(s)=x _(i)(s)(Ax(s))_(i)(x ^(T)(s)Ax(s))⁻¹.  (4)

The above replicator dynamics make sure that x₀≧0,x₁≧0, and x₀+x₁=1, andconstantly increases the average fitness value until it reaches thelocal maximum where each label becomes either 0 or 1. The diffusionprocess 103 encourages spatial homogeneity and adds disturbance to breaka tie in the reaction phase 104. After a few iterations, the outcome ofthe competition becomes clear and thresholding can be applied 106 at 0.5on x₁ to obtain the segmentation result 107. In experiments thediffusion and the reaction phases were applied alternatively for fouriterations 105. Other numbers of iterations may be implemented.

According to an embodiment of the present disclosure, a Gaussianellipsoid fit (EF) system implements a method which, given a markerpositioned near a target nodule, fits a 3D anisotropic Gaussian functionto the nodule's intensity distribution in a multiscale fashion. Anellipsoid that approximates the nodule's boundary is derived as aspecific equal-probability contour of the fitted Gaussian. Variousnodule size features (e.g., maximum diameter, volume, sphericity) aredetermined analytically from radii of the ellipsoid. The multiscaleanalysis is given by considering a Gaussian scale space of the inputsub-volume with a set of discrete analysis scales in an increasingorder, performing a Gaussian model fitting at each analysis scale, andfind the most stable estimate among the multiscale estimates byminimizing a form of Jensen Shannon divergence criterion. At each scale,Gaussian mean as the nodule center location is estimated by theconvergence of the scale mean shift procedures. In the neighborhood ofthe estimated mean, local data analysis is performed by the mean shiftprocedures initialized at a set of neighboring points. Gaussiancovariance matrix as the anisotropic spread is estimated by theconstrained least-squares solution of a linear system constructed withthe convergent mean shift vectors. Parameter settings may be selectedand varied depending on an application. For example, for the experimentsdescribed herein, the scale space is conceived with 18 analysis scaleswith 0.25 interval (0.5², . . . , 4.75²). And 35% 3D confidenceellipsoid is used for deriving an equal-probability contour from thefitted Gaussian.

The local distribution maximum applies thresholding at multiple levelsfollowed by connected component analysis on each thresholded volume. Itthen searches for object and their plateaus in the multiple thresholdedvolumes sequentially, starting from the one with the highest thresholdvalue in the order of descending threshold values. A new object is foundwhen a connected component has no overlaps to components in the previousvolume. An object becomes a plateau when the ratio between the volume ofthe object and the volume of its bounding box suddenly decreases by morethan some fraction (e.g., η) or the object merges with another plateau.Parameters of the method include the threshold values and η. η was setto 1/30 and the threshold levels were set to 0, 100, 200, . . . , 1100.Other parameter values may be selected.

The I-ELCAP protocol may be implemented for diameter estimation, whereinthe diameter of a nodule is estimated as the average length and widthwhere length is measured on a single CT image that shows the maximumlength from among a plurality of segments of a segmentation; and widthis defined as the longest perpendicular to the length of a segmenthaving the maximum length. The CT image is selected along the axialdirection, due to high-resolution and isotropic nature of the axialview. The I-ELCAP protocol may be implemented to estimate the diameterfrom the segmentation. Other measures of diameter may be implemented.

For LDM and RD, the diameter of a nodule is estimated as follows. Afterthe segmentation and connected component analysis, the component that isclosest to the marker is selected as the one corresponding to the noduleat the marker. In most cases, the marker is contained within thecomponent, but in some cases, it points to the background due toinaccuracy in marker positioning. Next, the component is analyzed sliceby slice in the axial view. For each axial slice, 2D connected componentanalysis is performed, an ellipse is fitted to each component, and thegeometrical mean of the axes is recorded for each ellipse. Among all 2Dconnected components, we select the one with the maximum geometricalmean for diameter measurement. The diameter is then estimated as theaverage width and height of the bounding box enclosing the selectedcomponent. Use of the geometrical mean instead of the arithmetic meangives a slightly better agreement with the manual measurements.

For EF, the diameter is estimated as follows. First, an ellipsoiddirectly derived from the estimated covariance is projected on the axialplane, and the radii of the ellipse on the projection plane isdetermined. The diameter can be estimated, for example, as thearithmetic mean of the radii times the scaling constant of 1.6416corresponding to the 35% confidence limit.

Referring to FIG. 1B, according to an embodiment of the presentdisclosure, for a hybrid method, when a nodule is attached to anotherstructure, RD segments both the nodule and structure together, resultingin a large volume. This segmentation volume usually stretches out to theboundary of the bounding volume (21×21×21 in our experiments). Thisobservation leads to a test for the nodule type: if the segmentationvolume touches one of six boundaries of the bounding volume, it isconsidered non-solitary. Otherwise, it is considered solitary. The ratiobetween the segmentation volume and the volume of the bounding boxenclosing the segmentation is determined to check if the segmentationhas a reasonably spherical shape. In experiments, only the boundarycheck was used to determine if the nodule is solitary.

By using this solitary/non-solitary check of a nodule with RDsegmentation, a hybrid approach to the diameter estimation problem canbe implemented. For each given marker of a VOI 111, a sub-volume of21×21×21 voxels were extracted, and RD segmentation 112 is applied tothe volume. A boundary check 113 is applied to the segmentation 112 ifthe boundary check indicates the nodule to be solitary, a diameterestimation 115 is applied on the RD segmentation 112. Otherwise, EF 114is applied on the sub-volume and the diameter is estimated 115. Thishybrid approach is denoted HB.

Referring now to the experiments discussed herein; the performance of HBhas been evaluated in estimating the diameter of pulmonary nodules. HBhas been compared against EF, LDM and RD. Centered at each marker placedby radiologists near a pulmonary nodule, a 21×21×21 bounding volume wasextracted and used as an input to the segmentation processes.

EF has difficulties in processing small nodules as the sample size forthe estimate is small. EF can also mistakenly include surroundingbackground for the estimate, which leads to overestimate of thediameter.

A problem associated with LDM is its sensitivity to the pre-determinedthreshold levels, which are typically set by a fixed increment. Thesegmentation of an object becomes inaccurate when the intensitydistribution of the object has an overlap with the distribution of itsplateau. The degree and the frequency of the problem can be reduced byincreasing the number of threshold levels, but at the cost of increasingthe computational load.

RD is effective in estimating the diameter of solitary modules.

The accuracy of the estimates were evaluated quantitatively by comparingthem with manual measurements by human experts (both radiologists andscientists in the medical imaging field). We use 1349 nodules taken fromover 240 CT volumes for the evaluation. The 1349 nodules are dividedinto 614 solitary and 735 non-solitary ones using the boundary test withthe RD segmentation. The performance of EF, LDM, and RD were evaluatedon solitary and non-solitary nodules separately. The accuracy ismeasured by the squared difference between estimates and thecorresponding manual measurements. FIGS. 2 and 3 show the results forsolitary and non-solitary nodules, respectively. In the plot, the meanand standard deviation of the squared difference errors are evaluatedamong nodules whose manual measurements are with ±0.5 of the centerdiameter (placed at a 0.5 mm increment between 1.0 and 7.0). The mean isshown in log scale to improve the visibility of the data while thedeviation is shown in linear scale. The purpose is to analyze theperformance in terms of the nodule diameter.

For solitary nodules, the error tends to increase with the centerdiameter of the analysis. This, it is more informative to compare theresults with the error normalized by the center diameter. RD is moreaccurate than EF and LDM across all sizes except at 3 mm where LDM with0.2913 normalized error is slightly better than RD with 0.2972normalized error. For non-solitary nodules, RD is the least accurate andEF is the most accurate one. For EF and RD, the error tends to decreasewith the center diameter, while for LDM, it stays around 2.5 mm². Theerror on RD clearly comes from segmentation of both nodule and anattached structure.

The performance of HB was evaluated using all 1349 nodules. FIG. 4 showsthe mean and the standard deviation of the squared difference error forEF, LDM and HB collected within the same diameter ranges as in FIGS. 2and 3. Both mean and deviation are shown in linear scale. HB constantlygives a smaller error than other methods in all sizes. This is predictedas HB uses RD and EF for solitary and non-solitary nodules,respectively, and the accuracy of RD and EF in their respectivecategories has been verified. Table 1 shows the mean and standarddeviation of the computation time for each estimation technique. Asshown, RD is computationally most efficient while EF is computationallymost expensive. HB is about 50% faster than EF. TABLE 1 ComputationalTime (seconds) METHOD EF LDM RD HB Mean 1.26 0.269 0.0510 0.651 Std.Dev. 0.25 0.18 0.014 0.7

It is to be understood that the present invention may be implemented invarious forms of hardware, software, firmware, special purposeprocessors, or a combination thereof. In one embodiment, the presentinvention may be implemented in software as an application programtangibly embodied on a program storage device. The application programmay be uploaded to, and executed by, a machine comprising any suitablearchitecture.

Referring to FIG. 5, according to an embodiment of the presentdisclosure, a computer system 501 for implementing a method forestimating solitary pulmonary nodule diameter with reactive-diffusionsegmentation can comprise, inter alia, a central processing unit (CPU)502, a memory 503 and an input/output (I/O) interface 504. The computersystem 501 is generally coupled through the I/O interface 504 to adisplay 505 and various input devices 506 such as a mouse and keyboard.The support circuits can include circuits such as cache, power supplies,clock circuits, and a communications bus. The memory 503 can includerandom access memory (RAM), read only memory (ROM), disk drive, tapedrive, etc., or a combination thereof. The present invention can beimplemented as a routine 507 that is stored in memory 503 and executedby the CPU 502 to process the signal from the signal source 508. Assuch, the computer system 501 is a general-purpose computer system thatbecomes a specific purpose computer system when executing the routine507 of the present invention.

The computer platform 501 also includes an operating system andmicroinstruction code. The various processes and functions describedherein may either be part of the microinstruction code or part of theapplication program (or a combination thereof), which is executed viathe operating system. In addition, various other peripheral devices maybe connected to the computer platform such as an additional data storagedevice and a printing device.

It is to be further understood that, because some of the constituentsystem components and method steps depicted in the accompanying figuresmay be implemented in software, the actual connections between thesystem components (or the process steps) may differ depending upon themanner in which the present invention is programmed. Given the teachingsof the present disclosure provided herein, one of ordinary skill in therelated art will be able to contemplate these and similarimplementations or configurations.

Having described embodiments for a system and method for estimatingsolitary pulmonary nodule diameter with reactive-diffusion segmentation,it is noted that modifications and variations can be made by personsskilled in the art in light of the above teachings. It is therefore tobe understood that changes may be made in embodiments of the presentdisclosure that are within the scope and spirit thereof.

1. A computer-implemented reactive-diffusion method for estimating a diameter of an object of interest comprising: providing a volume of interest comprising a plurality of voxels and including the object of interest; initializing at least two volumes of the volume of interest, wherein each of the plurality of voxels has at least two values corresponding to the at least two volumes of the volume of interest, respectively; performing, iteratively, a diffusion operation and a reaction operation on the plurality of voxels to adjust the at least two values; comparing, for each voxel, the at least two values to a threshold to assign each voxel to one of the at least two volumes, wherein the assignment of the plurality of voxels is a segmentation result; and estimating a diameter of the object of interest from the segmentation result, wherein the object of interest is represented by at least one of the at least two volumes but less than all the volumes.
 2. The computer-implemented reactive-diffusion method of claim 1, wherein the at least two volumes include a foreground volume and a background volume.
 3. The computer-implemented reactive-diffusion method of claim 2, wherein the at least two values correspond to a foreground label and a background label, respectively.
 4. The computer-implemented reactive-diffusion method of claim 1, wherein a number of iterations of the performing of the diffusion and reaction operations is predetermined.
 5. The computer-implemented reactive-diffusion method of claim 1, wherein performing the diffusion operation comprises performing a linear diffusion.
 6. The computer-implemented reactive-diffusion method of claim 5, wherein the linear diffusion is given by x _(i) Ý(s)=λ∇²x_(i)(s), where iε{0,1} and the diffusion rate A is set to 0.5.
 7. The computer-implemented reactive-diffusion method of claim 1, wherein performing the reaction operation comprises performing a competition between the at least two values in a form of replicator dynamics implementing a fitness matrix, wherein the replicator dynamics updates the at least two values.
 8. The computer-implemented reactive-diffusion method of claim 7, wherein performing the reaction operation comprises enforcing that each of the at least two values is greater than or equal to zero and the product of the at least two values is one.
 9. The computer-implemented reactive-diffusion method of claim 8, wherein the diffusion operation adds disturbance to break a tie between the at least two values in the reaction operation.
 10. A program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for a reactive-diffusion method for estimating a diameter of an object of interest, the method steps comprising: providing a volume of interest comprising a plurality of voxels and including the object of interest; initializing at least two volumes of the volume of interest, wherein each of the plurality of voxels has at least two values corresponding to the at least two volumes of the volume of interest, respectively; performing, iteratively, a diffusion operation and a reaction operation on the plurality of voxels to adjust the at least two values; comparing, for each voxel, the at least two values to a threshold to assign each voxel to one of the at least two volumes, wherein the assignment of the plurality of voxels is a segmentation result; and estimating a diameter of the object of interest from the segmentation result, wherein the object of interest is represented by at least one of the at least two volumes but less than all the volumes.
 11. A computer-implemented reactive-diffusion method for estimating a diameter of an object of interest comprising: providing a volume of interest comprising a plurality of voxels and including the object of interest; initializing at least a foreground volume and a background volume of the volume of interest, wherein each of the plurality of voxels has a plurality of values corresponding to the foreground volume and the background volume, respectively, wherein the plurality of values indicate whether a voxel is likely to belong to the foreground volume or the background volume; performing a diffusion operation and a reaction operation on the plurality of voxels; updating the plurality of values; comparing, for each voxel, the plurality of values to a threshold to assign each voxel to one of the foreground volume or the background volume, wherein the assignment of the plurality of voxels is a segmentation result; and estimating a diameter of the object of interest from the segmentation result, wherein the object of interest is represented by the foreground volume.
 12. The computer-implemented reactive-diffusion method of claim 11, wherein the diffusion and reaction operations are performed iteratively for a is predetermined number of iterations.
 13. The computer-implemented reactive-diffusion method of claim 11, wherein performing the diffusion operation comprises performing a linear diffusion.
 14. The computer-implemented reactive-diffusion method of claim 1, wherein performing the reaction operation comprises performing a competition between the plurality of values in a form of replicator dynamics implementing a fitness matrix, wherein the replicator dynamics updates the plurality of values.
 15. The computer-implemented reactive-diffusion method of claim 14, wherein the diffusion operation adds disturbance to break a tie between the at least two values in the reaction operation.
 16. The computer-implemented reactive-diffusion method of claim 11, wherein each of the plurality of values is greater than or equal to zero and the product of the plurality of values is one. 